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ABSTRACT 


A dynamic structural model of a submerged ring is 
developed using trigonometric series. It is constructed 
for use in cOnjunction with a finite element fluid model 
to examine the effects of cavitation on underwater shock 
loading of a structure. The governing equations and the 
time integration algorithm used in the model are described. 
Results predicted by the model are compared to known 
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results. The program listing 1s given. 
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I. INTRODUCTION 


The program evolved during this study models a submerged 
Circular ring using trigonometric series. It was developed 
for use in cOnjunction with a finite element fluid model 
based on a displacement potential formulation. The purpose 
of the combined models is to predict the effects of 
cavitation on underwater shock loading of the structure. 
More information on the fluid model can be found in 
References 1 and 2. The purpose of this paper is to 
describe the development of the structural model and to 
present some of the results obtained from its use in 
combination with the fluid model. A listing of the FORTRAN IV 
program implementing the structural model 1S given in the 


Appendix. 
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II. GOVERNING EQUATIONS 
POL PPERENTIAL EQUATIONS FOR THE DEFLECTION 
OF A THIN CIRCULAR RING 
In Figure 1 the solid line represents the ring after 
' deformation, and the dotted line the ring before deformation. 
For small deflections the curvature of the element m,n, can 


be taken as 


ag _ dé@ + Adé 
Ry ds + Ads 


where R, w, 9 and s are defined as shown in Figure 1 and 


w 1s taken positive inward. Making use of the relations 


2 2 
Ado = de, - de = 2+ S¥a-S = Sas 
ds ds 
ds 
Las = ds, - ds = (r-w)d@ ~ rd6 = -wd? = --wW z 
and substituting into the equation above yields 
De 2 
daw dw W 
a ds _ R + ds 
l W _ W W W W 
ds (1 - 5) ds (1-5) (1+) ds (1 m) (1+) 


Neglecting the higher order terms, the above expression 


becomes 


he 
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Figure 1. Geometry of Ring Deformation 
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aD 2 
al dé W dw il W déw 
— «= = (l +s) + — = =(Ll+ ez) + 
Ri ds R as R R ds? 
or alternately 
LE 2 

—- = == + ex = > (w + 2) (1) 
ik R ds R dé 


For a ring where the thickness is small compared to the 
radius and elastic behavior is assumed, it can be shown 
that the approximate relationship between deflection and 


loading is [Reference 7] 


) 
tt 
{ 
OS 
Nw 


ally 
So 


where M is the bending moment about the centroidal axis 
meets the flexural rigidity of the ring. A positive 
bending moment produces compression in the outside fibers 
of the ring. Combining equations 1 and 2 yields the 
differential equations for the deflection of the ring given 
below. 

Oe. yj a 


R 36° 


(is 
w 


14 





B. GOVERNING EQUATIONS OF MOTION 
The governing equations of motion used in the program 

were arrived at by the application of Hamilton's principle. 
iemergder tO apply Hamilton's principle, it is first necessary 
to derive expressions for the strain energy and kinetic 
energy of the ring as well as the work done by the external 
loads. In these derivations the ring was taken to be 
homogeneous, elastic, and of unit width. The pressures on 
the ring and its deflection were represented by the 
pressures at and the deflection of a set of nodal points 
equally spaced along the circumference of the ring. 

. The shock front is assumed to approach the ring normal 
to the 9 = 0 plane as shown in Figure 2, where 9 is taken 


to be positive counterclockwise. 


Ligeia mince ohnock Front Orlentation 


ie 





———————— ————— — 


Since the ring iS symmetric about the 6 = 0 plane, 
results are given for 9 from 0° to 180°. The response of 
the ring is represented by the response of a set of N+l 
nodal points equally spaced along the ring. The radial 
displacement 1S approximated by a trigonometric cosine 


series 
a eeeceecos oc). (4) 


The tangential deflection is represented by the trigonometric 


Sine series 
i een cl (5) 


where v is taken to be positive in the negative 8@ direction. 
Similarly the normal pressure applied to the ring is 


represented by 


a = » Cc, cos né . (6) 


ime Strain Energy 


The strain energy is comprised of two components 
[Ref. 4]. One component is the strain energy due to bending 
and the other to strain in the 8 direction (E,)- The strain 


energy stored in the ring due to bending is calculated as 


16 





the work done by the internal moment acting over the 
emiermte ring, and for a ring of unit width is 


l 27 
U = z i M RKdé 


where k is the curvature of the ring. Combining this 
equation with equation 3 and using the relation «x = M/D 


yields 


2a 2 
3 i f py AB+w* ae (7) 
Ree 


The strain energy due to the strain in the §$ direction is 
Obtained as the distance traveled by the average normal 
force in the @ direction. 

The circumferential strain 1s composed of two 
components the first of which is the result of the radial 


displacement of the ring (=) and the second results from 
Jl veNt 


tangential displacement and is (= 5@) - The strain energy 
Seam the normal force (N) is 
l 27 l Png 
U, = 5 J Nek do = = f Ae,Ee,R dO 

0) 0) 
Sha (8) 

25 27 

_ en 2 eT AE dv 2 
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where A is the cross sectional area of the ring and 


Eis the modulus of elasticity of the ring material. 


Combining equations 7 and 8, the expression for the 


total strain energy is 


Substituting into this equation the series expressions for 
radial and tangential displacements, and assuming constant 
geometrical and material properties over the ring yields 


after integration 





An expression for the change in strain energy 6U 
corresponding to small changes in tangential and radial 


displacements (oR andséa_) Can now be found 


wD 


N 
su = -—= [ } me 1) 4 éa_] + al 
n=0 yh 


N 2 
5 = [ d (n be + na) ob, 
n=0 


- (nb _ + a) oa.) | (10) 


* 

In this equation and several that follow which contain 
a summation expression starting with n = 0, a factor of 2 
associated with the n = 0 term is omitted for the sake of 
brevity. The term is accounted for in the solution process 
used with the model. 
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2. Work 
The work associated with a small change in radial 
displacement 1S given by 
277 


jw = ‘i p ow R dé 
0 2 


Substituting the relations 


N 

p = ) Cn COS mo 

m=0 
and 

N 

ébw = ) 6a cos né 

n 

n=0 


Pn kinetic Energy 


The expression for the kinetic energy of the ring 
can be arrived at by considering an infinitesimal section 
of the ring that has been set in motion. The kinetic energy 
of the element can be represented as the sum of a contribu- 
tion due to translation of the mass center and a contribution 


due to rotation. This leads to 


9 





—_ Sele ——eE 


a ae ae a Iw, 2 
dT = 5 9(ARd6) (wot vi) + 5 eIRd6 (aaq) 


where I is the area moment of inertia about the centroidal 


axis. 
Using the series relations 
N N 
w = } a_ cos né6 a) acoso 
n=0 a n=0 n 
N N. 
oe) b, Sin ne a ne 
n=] n=1 
yields 
1 ss 2 a 2 
See > ORGS) I( ) a, cos no)” + ({ } b. sin né)"] 
5 n 
n=0 =] 
1 _ 2 
1 Ge @abistells) (3 ae Sane G ) 
n=l 
Integration over the ring yields 
7 uss) y i (UR eae A tls er 
2 a6 R n 


The resultant change in kinetic energy due to a 


small change in velocity components can then be found as 


sT = mo } [(AR + n* 2)a 6a, + ARb 6b] (12) 
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mero calion OL Ham ton’ s Principle To Find 
GCeupled Equations Of Motlon 


From Hamilton's principle [Ref. 5] it is known 


foac 


f (ie 60 + ow) idt = 0 jab} 


Semomiing equations 10, 11, 12 and 13 yields, after 


Carrying out an integration by parts on the first term, 








. Det: , a 
) mo(AR + n“ =) [ada - f a é6a_ at] 
& She > Gialaiieg) ial 
n=0 s ty 
1 
N i, a 
+ y TPAR[D sb, | = f b., Sb dt] 
Bo t a 
1 
\s N N 
2 
- f*@2 ey wP-1)* a, 6a) + SBE OF in? + nal) db 
1 r= 
“ 
gmnb. + a )éa_ lJ = } c. éa mR} dt = 0 
n n n & n n 
n=0 
Since the 6a, and 6b, can be chosen Sreokenarilvyey, this 
becomes 
eo. c, 2 r,2, 2 2 BO es ah 
ie + (=) eee le =) ete = 1) ees os) a an 
Cac Cee. 2 a 
be + (3) na, + (3) n oo = 0 (14) 
ce 2s | 
where c = (E/p) is the speed of sound in the ring and 


ZA: 





where r is the centroidal radius of gyration of the ring 


section and the relation I = Ar? has been used. 


See ocparated Equations of Motion 


It is common practice to separate ring deformations 
into flexural modes and extensional modes. The initial 
version of this program utilized such a resolution. The 


flexural mode is characterized by the requirement that 


This implies that the Fourier coefficients satisfy the 


relation 


a + nb = 0 Sy), 


The extensional mode is defined to be geometrically orthogonal 


to the flexural mode so that 


na - b = QO (16) 
(E) 





2 
[oA( (= =) + zn )ja + (=5 aha = 3 a, 7a 
n (EF) R ic) 
b = -a foe (a7) 
(FP)  (F) 


oe 





If equation 16 is substituted into equations 14, the result 


is 


Initial solutions were carried out using equations 


17 and 18 and then combining the results by 


This procedure gave results which satisfactorily predicted 
the ring bending moments but gave poor accuracy for axial 
force determinations. 

Re-examination of the foregoing solution process 
revealed that the mode shapes defined by equations 15 and 
16 are not orthogonal with respect to the mass or the 
stiffness matrices of the system. It was accordingly 


decided to return to equations 14 and solve them 
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Simultaneously in the time integration algorithm. This 
change led to accurate axial force determinations. The 
separated equations of motion, 17 and 18, are used in the 


program as described in Section III, b. 
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Ili. TIME INTEGRATION METHOD 


A. CENTRAL DIFFERENCE ALGORITHM 
Time integration is accomplished using a central 


difference algorithm. If ans bu and Cc, are known at 


time p fh), then at) (the value of a at t {2), may be 


evaluated from the first of equations 14. Letting h 


(1-%3) 


represent the time step and a represent a at time 


t 2) - h/2, the next value of a 1s calculated from 


Aer : pV ee 
“led = 4° a) 4 na \?? (19) 
Using this result, the next a, (at time fo = et) + h) 
memcalculated from 
Bute 
4 (ttl) — 4‘) a h peitea) (20) 


n pl gl 


The value Dh 1S similarly found from the second of 
equations 14. A pair of equations paralleling 19 and 20 


| Bo ep i 
are used to calculate ia ese 


and b When these steps 
have been completed for each n, the value of radial 
displacement w is found at each structural node from 
equation 4. These displacements are passed to the fluid 
program where they furnish required interface boundary 


(1+1)- Values of 


conditions to allow an advance to time t 
fluid pressure at the interface nodes are returned by the 


Guid program. 


2 








The new nodal pressures are used to calculate c,'s at 
ae). It is then again possible to advance an and a. 
using equations 14, 19, and 20 and to perform the parallel 
calculations for b. and bo: 

The solution process is started with the ring at rest 


under uniform hydrostatic pressure. Under the loading ay 


is the only nonzero Fourier coefficient. Because at = 0 
eee te 
is known (rather than at a) a fictitious starting value of 
Si be 
at SMe first calculated from 
ese - 
gi 5) ae Ah 3 (0) 
n 59) 


A similar starting procedure is used for bs. 
Een LECTION OF THE NUMBER OF 

VIBRATIONAL MODES USED 

The number of vibrational modes that can be modeled 
accurately is limited by the numerical integration algorithm. 
Specifically, the algorithm becomes unstable for time steps 
in excess of about 0.3 of the period of the structural mode. 
‘The accuracy of the algorithm deteriorates even before the 
stability limit is reached. For this reason, a criterion 
was established for limiting the number of modes used, 
based on the time steps selected. 

In the program the separated equations of motion for 
the extensional and flexural modes were used to find the 


frequencies needed in applying the criterion. From 
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equations 17 and 18 the natural frequencies of the 


extensional and flexural modes were found to be 


wi = (€)* (1 + n°) 2 ee 
(E) 
2 2 
vw? = (S) * (Ey n es =i) ae 
(F) mel 


The criterion used was that the period for the highest 
mode retained be at least five times the time step used. 
Since the extensional modes have higher frequencies for a 
mewn, ENLS Criterion comes into effect first for the 
extensional modes. Once the point is reached (where 
n = n, (max) ), the program switches from the coupled equations 
of motion 14, to the separated equation 17. When this 
occurs the higher extensional mode coefficients (n > n, (max) ) 
are not needed, and only the flexural mode coefficients 
for n > n, (max) ) are used. If the time step is large 
enough that the criterion is also met by the flexural 
modes (n = n,, (max) ) then both the a's and b_'s are omitted 


for all modes where n is greater than n,, (Max) . 


AT 





ye neo UELs 


To check the accuracy of the results obtained from the 
program, two kinds of comparisons were made. First, the 
solution of equations 14 arrived at by the program were 
compared to an analytic solution of equations 14. Second, 
to verify that equations 14 correctly predict ring behavior, 
calculated dynamic response of the ring to two particular 
loadings was compared to known static results [Reference 6] 


for the same loadings. 


A. COMPARISON OF MODEL RESULTS TO EXACT SOLUTION 
The check on solution accuracy was based on the exact 
solution for a ring initially at rest and suddenly loaded 


by a steady pressure 


P = [1.625 - .25 cos 8 - .375 cos 28] MPa 


The exact solution, using the parameters specified in 


Table IT, is 
ay = Aerie om tit eneos HC1O1] . 203%) | 
ay = Go, ae ~ LBs ieelo im - cos (1429.t)] 
a. = -.104129[{1 - cos (85.616)t] 
5 


- 3.7355xl0 ~{l1 - cos (2259.7t)] (22) 
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Table I gives results for radial displacements at 
t = 1 msec for 6 = 0°, 90° and 180°. Computer results are 
given for four different time steps. Even the coarsest 
time step gives results within 1% of the exact values 


Obtained from equations 22. 


B. COMPARISON OF MODEL RESULTS TO STATIC BEHAVIOR 

To compare the dynamic response predicted by the program 
to the known static response of the ring for the same 
loading, two simple cases of loading were assumed. The 
first loading case examined was uniform pressure surrounding 
the ring; the second case was two equal and opposite concen- 
meeeea JOads acting 180° apart as shown in Figure 3. The 
second loading case could not be represented exactly by 
the finite trigonometric series loading representation used 


by the program. It was approximated by nonzero values of 


Seeesure Only at 6 = 0° and @ = 180°. 


Figure 3. Concentrated Loading on Ring 
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TABS Eaae 







mxl0- mxl0- mxl0- 


ee 

— | pe 

eee f= p= 
equations 17 


Ring Parameters: a 5m, A= gee m* , B= 200 GPa 
= 7830 kg/m?, r = .158m 





Comparison of program solution with 
analytic solution 


TABLE it 






Twice Known -10.0 
Static Results 


Ring Parameters: 5m, A= ue m? , E= 200 GPa 
250 kg/m fee Seam 





Comparison of program dynamic solution 
with twice known static solution for 
uniform external pressure on ring 
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Comparison between the dynamic and static responses 
of the ring were made for the radial displacements, bending 
moments, and forces at particular nodal points on the ring. 
For the uniform pressure case any nodal point on the ring 
can be used since the displacements are uniform over the 
ring. For the concentrated loading case, the nodal points 
chosen were at the points of application of the forces and 
midway between the application points. 

For the dynamic response of the ring predicted by the 
program, the ring is taken to be initially at rest and 
undeformed. The anticipated behavior of the ring after 
the application of the loading is for the ring to deflect 
through the radial displacements found for the static 
loading and to continue to deflect until a radial displace- 
ment with a magnitude of approximately twice the static 
displacement is reached. The maximum deflections were 
expected to occur after a time approximately equal to 
one half the period of the vibrational mode that dominates 
the deflection of the ring. For the uniform loading case 
the dominant vibrational mode is the fundamental extensional 
mode. For the concentrated loading case, it is the second 
Flexural mode. The periods for these vibrational modes can 
be found using equations 21. Tables II and III give the 
comparisons between the dynamic response from the program 
and twice the known static response for the two loading 


cases considered. As is seen, relatively good correlations 
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TABLE III 


(max 
i ae Hea 
at B 


Model 
-2965 ~2839 3.760 P2910 


Static 

Analysis a = 
ee .2634 a7). 85536 
Times 2 


Model 
1 | zs] -.t887 | -2. 383 
jdysts 2 ee LS So. 28/0) -4.543 


Model 
e510 = eles OS Sa aS -3.544 


~ 1425 rege lemre NES: oe) 


Sma = 05-7, = = 200 GPa 
7830 kg/m3, aca FIs 





Ring Parameters: 


il ‘i 


Comparison of dynamic solution with twice known 
StamlewscGIntlOn ftOm~elOaging Of Fig. 3 


BZ 





between the dynamic responses and twice the static responses 
were obtained. The maximum values shown in Table III occur 
at approximately 35.3 milliseconds after time zero, which 
compares favorably to one half of the period of the second 
flexural mode (36.6 milliseconds). The maximum values in 
Table II occur at 3 milliseconds which is approximately 

one half the period of the fundamental extensional mode 


(3.11 milliseconds). 


Cs OUTPUT EXAMPLES 

A sample of the output produced by the structural 
program is shown in Figures 4 and 5. The print code at the 
top of Figure 4 is an input code used to select the informa- 
tion desired in the print out. Following the print code 
are two sets of input parameters to be used in the structural 
model. These are cotlowed by the mass and stiffness 
coefficients calculated for each of the vibrational modes. 
This information remains constant for any given run and 
1s therefore only printed once. The rest of the information 
varles with time and can be printed as often as desired. 
The pressures at nodal points generated by the fluid model 
are printed out; these pressures are for one time step 
before the time printed below the pressures. Following the 
time come the Fourier coefficients for the radial displace- 
ments, tangential displacements, and pressures. In 


Figure 5 at the top, the bending moment and the axial force 


3.3 





at each nodal point are printed for the time given in 
Figure 4. Tabulation of the radial displacement follows. 
Figure 6 is an example of the graphical output obtained 
From the fluid model when the fluid and structural models 
are run jointly. It is a time sequence Ge elght  oloOts Of 
Fluid nodal pressures over the domain. The plots are shown 
at 8 ms intervals. The left hand edge of each plot repre- 
sents the plane of symmetry and the top and bottom rows 
represent, respectively, the entry and exit faces for the 
shock. The x's on the lower left side of the plots 
represent dummy nodes inside the structure. The pressure 
ranges for the mapping characters are shown in Figure 7. 
In Figure 6 the development of the cavity cantle followed. 
It develops in the second, third, and fourth frames, 


collapses in the fifth, and has vanished by the sixth. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


The comparisons made between the program results where 
convergence of the solution has been obtained to results 
obtained by other methods, show that the model satisfactorily 
predicts the dynamic response of a ring to external pressures. 
The ring model cannot predict exactly the behavior of a real 
structure such as a submarine hull. However, by a suitable 
choice of ring parameters the response of the ring structure 
will provide a good first approximation to the response of 
a real structure. 

A useful improvement to the program would be to include 
a section which would automatically pick out maximum values 
of bending moments and axial forces in the ring as well as 
when and where they occur. This would allow any desired 
time intervals between data printout to be used without 


missing peak responses. 
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